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ABSTRACT 

A  numerical  procedure  Is  presented  for  the  solution  of  the 

characteristic  equations  for  a  non-lsoenergetlc,  supersonic,  real  gas 

jet  expanding  Into  a  constant  pressure  ambience.  Included  Is  a  procedure 

for  Inserting  the  expansion  fan  and  the  embedded  shock.  Examples  are 
given  for  high  temperature  argon  jets  with  radiation. 
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NOMENCLATURE 
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nr  Number  density  of  the  rth  specie 

p  Pressure 

0  Flux  of  dissipating  energy 
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Rotational  partition  function 
Translational  partition  function 
Vibrational  partition  function 
Radial  coordinate 

Transformed  coordinate  tangent  to  oblique  shock 
Upstream  shock  point 
Downstream  shock  point 
Temperature 

Axial  component  of  velocity 
Volume 

Virtual  point 

Radial  component  of  velocity 
Magnitude  of  velocity  vector 

Right  hand  side  of  various  finite  difference  equations  as 
defined  In  text 
Electronic  partition  function 
Ionic  charge 
Axial  coordinate 
Oblique  shock  angle 
Density 

Local  flow  angle 

Characteristic  temperature  of  rth  species 
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CHAPTER  I 
INTRODUCTION 

I.  PURPOSE 

With  the  increased  research  emphasis  on  high  temperature  jets, 
the  effects  of  radiation  as  well  as  real  gas  effects  have  taken  on 
Increased  Importance.  Consequently  the  analytical  tools,  l.e.  the 
mathematical  models,  for  studying  these  jets  need  further  development. 

Available  solutions  make  at  least  one  If  not  all  of  the 
following  assumptions:  (I)  the  gas  is  perfect;  (2)  stagnation  enthalpy 
Is  constant  along  streamlines;  (3)  the  flow  is  I rrotational ;  (4)  embedded 
shocks  can  be  neglected. 

The  purpose  of  this  study,  then,  is  to  solve  the  system  of 
characteristic  equations  obtained  without  any  of  these  assumptions. 

It  is  anticipated  that  the  resulting  computer  program  will  be 
a  valuable  research  tool. 

II.  PHYSICAL  CONSIDERATIONS 

When  a  supersonic  jet  exhausts  into  a  lower  pressure  environ¬ 
ment,  a  complex  flow  field  is  created  as  In  Figure  I. 

The  fluid  at  the  nozzle  lip  goes  through  an  immediate  expansion 
to  the  ambient  pressure.  Since  In  supersonic  flow  disturbances  are  prop¬ 
agated  downstream  only,  the  remainder  of  the  fluid  in  the  same  plane  as 


1 


AEDC-TR-69-1 84 


the  nozzle  exit  remains  undisturbed.  However,  as  the  core  of  this  jet 
encounters  the  rarefaction  propagated  from  the  nozzle  lip.  It  too  goes 
through  an  expansion.  When  the  flow  encounters  the  rarefaction  propagated 
across  the  centerline  from  a  point  on  the  nozzle  lip  diametrically  oppo¬ 
site,  further  expansion  Is  Induced  such  that  a  pressure  much  lower  than 
ambient  Is  reached.  The  flow  must  then  recompress  to  the  boundary  pres¬ 
sure.  Thus,  a  pressure  trough  could  be  said  to  form  that  deepens  as  the 
core  of  the  flow  encounters  the  overlapped  rarefactions.  While  near  the 
boundary,  the  pressure  remains  close  to  ambient.  Until  It  encounters 
the  higher  pressure  region  near  the  boundary,  the  core  Is  not  Influenced 
by  whether  It  Is  expanding  Into  a  modest  pressure  difference  or  Into  a 
vacuum. 

It  is  convenient  to  think  of  the  field  as  consisting  of  two 
regions:  an  upper  and  a  lower  field.  In  the  lower  field  occur  the 
combined  rarefactions  causing  the  overexpansion.  It  is  tempting  to  call 
this  a  source  flow  field  as  Indeed  It  is  assumed  in  Reference  (1)*. 

That  is,  the  flow  appears  to  radiate  from  a  point  source. 

At  odds  with  this  Is  the  upper  field  which  is  at  a  higher 
pressure  since  it  is  Influenced  by  disturbances  from  the  boundary. 

Since  the  boundary  of  the  flow  Is  a  streamline,  the  dominant  direction 
of  the  upper  field  could  be  considered  as  being  axial. 

Obviously  along  some  line  in  the  plane  being  considered  there 


^Numbers  in  parentheses  refer  to  similarly  numbered  references 
in  the  bibliography. 
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will  be  a  conflict.  On  this  line  multiple  conditions  appear  likely  tq 
occur.  The  flow  will  at  once  encounter  a  dominantly  source-1  Ike,  low 
pressure  region  and  a  mainly  axial,  higher  pressure  region. 

An  oblique  shock  serves  to  couple  the  two  differing  fields. 

This  shock,  variously  named  an  embedded  shock,  boundary  shock,  secondary 
shock,  barrel  shock,  and  incident  shock,  starts  as  a  weak  shock  emanating 
from  the  nozzle  lip  following  the  outer  portion  of  the  expansion  fan.  As 
the  pressure  trough  deepens,  the  shock  strengthens  and  finally  curves 
toward  the  axis.  The  method  of  characteristics  solution  predicts  the 
Intersection  of  the  shock  with  the  axis  If  subsonic  conditions  are  not 
encountered  first.  However,  in  actuality  the  phenomenon  of  the  Mach  disc 
occurs  before  then.  The  Mach  disc  Is  a  normal  shock  behind  which  there 
Is  subsonic  flow. 
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CHAPTER  II 
THEORY 

I.  DERIVATION  OF  CHARACTERISTIC  EQUATIONS 

The  conservation  equations  are  derived  in  most  textbooks  on 
theoretical  gas  dynamics.  For  an  axially  symmetric,  inviscid,  compress¬ 
ible,  non-isoenergetic  flow  they  are  considered  as: 
for  conservation  of  mass. 


=  0 


(I) 


for  conservation  of  axial  momentum. 


3u  3u  ,  3 p  n 

pu  37+pvTF  +  5i’-° 


3u 


(2) 


for  conservation  of  radial  momentum. 


pu  3v  3p 

3z  k  3r  3r 


=  0 


(3) 


and  for  the  conservation  of  energy. 


pu  +  pv  1^.  =  Q(z,r,u,v,p,h) 


(4) 


This  is  a  set  of  four  equations  for  the  four  unknowns  u,  v,  p,  and  h  since 


H  =  ^(u2  t  v2)  +  h 


(5) 


and  from  the  equation  of  state. 
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p  =  f(p,h)  (6) 

Following  Reference  (2),  Equations  I  to  4  are  rewritten  In  a  form  such 
that  only  the  dependent  variables  appear  In  the  derivatives  with  respect 
to  the  coordinate  variables.  They  become,  using  Equations  5  and  6, 


(7) 


3p  3h  3p  3h 

+  n»  air  97  +  ™  sk  W 


=  -  pv 


pu 


3u  ,  .  3u  .  3p  _  A 

57  +  pv5F+5f-  0 


(8) 


3v  3v  3p 

pugj  +  pv  7F  +  3r  =  0 


(9) 


pu2  lr+  puv  If+  puv  Ii+  pv2If  +  pu  §f  +  pv  jfs  q 


(10) 


These  four  equations  along  with  the  definitions: 


3u  .  .  3u  ^  . 

77  dz  +7r  dr  =  du 


(II) 


£**&*»* 


(12) 


|f  dz  +  dr  =  dp 


(13) 


ill.  dz  +  |ll  dr  =  dh 
3z  3r 


(14) 


are  considered  to  be  a  linear  system  of  eight  algebraic  equations  and 
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eight  unknowns.  The  unknowns  are  taken  to  be  the  partial  derivatives  of 
the  dependent  variables  with  respect  to  the  Independent  variables. 

The  question  Is  then  asked,  "Do  curves  exist  within  the  domain 
of  Interest  on  which  these  derivatives  may  be  discontinuous?" 

it  can  be  shown  that  these  curves  do  exist  if  the  function  Q 
In  Equation  10  is  not  a  function  of  the  derivatives  of  the  dependent 
variables. 

This  involves  a  tedious  expansion  of  large  determinants  which 
will  not  be  presented  here.  The  resulting  equations  are 

udr  -  vdz  =0  (15) 

which  Is  recognized  as  the  equation  of  a  streamline; 

(u2-a2)dr2  -  2uvdrdz  +  (v2-a2)dz2  =0  (16) 


where  because  of  the  choice  of  dependent  variables 


3p  .  .  3p 
W+  p  3^ 


and 


(udr-vdz) (udv-vdu)  +  (udr-vdz)z  v/r 


+  (vdr+udz)  t  (udr-vdz)2  Q  =  0 


which  holds  along  the  curves  described  by  Equation  16; 

also 

dp/p  +  udu  t  vdv  =  0 


(17) 


(18) 


(19) 
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which  is  Bernoul 1 1 *s  equation;  and 

dh  +  udu  +  vdv  =  Qdz/pu  (20) 

which  is  recognized  as  the  energy  equation  aiong  a  streamline. 

The  curves  described  by  Equation  16  are  referred  to  as  the 
left  or  right  running  characteristics,  while  Equation  17  is  simply  the 
acoustic  speed.  Equations  18  to  20  are  usually  referred  to  as  the 
compatibility  relations. 

From  the  original  conservation  equations,  a  rather  remarkable 
set  of  relations  has  been  obtained.  These  are  relations  for  three 
families  of  curves  along  each  of  which  is  given  the  relationship  among 
the  flow  properties  In  the  form  of  ordinary  differential  equations. 

These  equations  are  put  In  finite  difference  form  and  solved 
numerically.  The  detailed  procedure  is  presented  in  Chapter  III. 

II.  EXPANSION  AT  THE  NOZZLE  LIP 

At  the  nozzle  lip  the  flow  must  turn  abruptly  to  match  the 
ambient  pressure.  In  an  inviscid  flow,  the  flow  properties  are  not  single 
valued,  I.e.  the  field  at  that  point  Is  assumed  to  experience  a  con¬ 
tinuous  change  In  pressure  from  the  pressure  existing  at  the  nozzle 
exit  to  the  ambient  pressure.  A  relation,  then,  is  needed  for  the 
properties  at  a  point  as  a  function  of  pressure. 

This  relation,  the  so  called  Prandtl -Meyer  relation,  can  be 
derived  in  many  ways,  it  is  interesting  to  see  that  It  can  be  obtained 
readily  from  the  compatibility  relations. 
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Dividing  Equation  18  through  by  udr-vdz  and  taking  the  limit 
as  dr  and  dz  go  to  zero  along  a  left  running  characteristic,  the 
expression 


I  imlt 
dr  •+  0 
dz  •*  0 


r 


udv  -  vdu  t  (udr  -  vdz)  —  + 

r  p 


vdr 


udr 


t  udz  ] 
-  vdz  J 


becomes 


+  -U-  <udr  -  vdz)Q 

p2  9h 


udv  -  vdu  t 

P 


fllmit 
dr 
dz 


•»+  r  W 

■+■  0  vdr  +  udz  f  =  o 
0  [udr  -  vdz  l 


(21 ) 


By  taking  the  limit  along  the  path,  the  slope  of  which  Is  given  by  Equa¬ 
tion  16,  i.e. 


dr  uv  +  ay  u2  -  a2  +  v2 
dz  =  u2  -  a2 


(22) 


Equation  21  becomes 


udv  -  vdu  t 


=  0 


(23) 


where 

m2  _  u2  +  v2 

•-7- 


Also,  the  limit  of  the  compatibility  relations  along  a  streamline  is 
taken,  obtaining 
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and 


udu  +  vdv  = 


.  & 
P 


(24) 


dh  +  udu  +  vdv  =  0  (25) 

Equation  25  Is  Integrated  directly  giving 

2  2 

h  +  - — ^—V  =  constant  (26) 

Equations  23  and  24  are  then  arranged  such  that 


dv 

(Ip  “  v  (Ip  “ 


v 


du 

u  Ip  = 


Jl 

P 


Using  Cramer's  rule 


du 

dp 


u  /m2  -  | 
F+  v  ' '  'p  - 

u2  +  v2 


and 


dv 

3p 


2  2 
U*  +  V^ 


(27) 


(28) 


Equations  27  and  28  can  be  written  more  familiarly  In  polar  form  for  the 
velocity  as 
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de  7m2  -  i 
■ap  =  “  &r- 

and 


(29) 


dp  pw 

Thus,  wl+h  Equations  26,  27,  and  28,  the  variables  u,  v,  and  h  can  be 
determined  as  a  function  of  p.  Consequently,  the  properties  at  a  point 
on  the  nozzle  lip  are  determined  by  varying  pressure  between  Its  Initial 
value  and  the  boundary  value  and  solving  for  u,  v,  and  h. 


III.  EMBEDDED  SHOCK 

As  discussed  In  the  Introduction,  an  embedded  or  secondary 
shock  Is  encountered  In  the  flow  field.  In  the  physical  field  the  shock 
occurs  when  the  Mach  waves  coalesce.  The  shock  wave  occurs  In  the 
solution  field  when  characteristic  curves  of  the  same  sense,  say  right¬ 
running,  coalesce  and  finally  overlap  one  another.  This  Is  treated  as  a 
discontinuity  In  the  dependent  variables.  To  cross  the  discontinuity  the 
oblique  shock  relations  must  be  utilized. 

The  shock  relations  fellow  from  the  conservation  equations. 
Equations  I,  2,  3,  and  4.  Consider,  for  example,  the  energy  equation. 
Multiplying  the  continuity  equation  by  H  and  adding  it  to  the  energy 
equation  gives 


3(puH)  .  3(pvH)  „  pvH 
3z  3r  =  9  "  r 
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Now,  consider  the  surface  of  discontinuity  of  the  dependent 
variables  as  occurring  at  an  angle  B,  Figure  2,  with  the  axial  coordinate 
at  the  point  of  interest  In  the  r-z  plane.  Transforming  Equation  30  In 
the  r-z  plane  to  coordinate  normal  and  tangent  to  the  curve,  say  n  and 
s.  Equation  30  then  becomes 

sin  a  +  3(^H)  cos  a  +  cos  a  _  3(^H)  sin  a  =  q 

_  pvH  (31) 

r 

Adding  terms  to  both  sides  to  complete  the  derivative  gives 


•  (puH  sin  a  +  pvH  cosa)  +  J_  (puH  cos  a  -  puH  sin  B> 
3n  3s 


=  0-4^-  <PvH  cos  ®  +  s,n  B)  - 

ds  r 


(32) 


Taking  the  integral  of  both  sides  with  respect  to  n  gives 


i: 


Nr 


(puH  sin  a  +  pvH  cos 


a)  dn  +  /  J_  (puH  cos  a  -  pvH  sin  a)dn 
m  9S 


72 


Qdn 


72  "2 

-J  (pvH  cos  a  +  puH  sin  a)  —  dn  -  I 
m  ds  -V,  r 


dn 


(33) 


Assuming  the  existence  of  the  Integrals  and  the  continuity  of  the 
Integrand  of  the  first  Integral  on  the  left  hand  side  for  ail  N(  -  n  - 
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(puH  sin  g  +  pvH  cos  8)  ^  -  (puH  sin  6  +  pvH  cos  8)|nj 


r 

+  /  I- 

3s 


(puH  cos  8  -  pvH  sin  8>dn  = 


Qdn 


-r 


(pvH  cos  8  +  puH  sin  8)  —  dn  -  . 

ds  J 


pvH 


dn 


(34) 


I 


The  Integrals  then  vanish  In  the  limit  as  N2  and  N|  are  made  to  approach 
a  point  on  the  discontinuity  from  opposite  sides,  or 

P2u2H2  s*n  B  +  P2V2H2  cos  ®  “  plul^l  s*n  &  +  PIVIH|  cos  ®  (35) 

A  similar  procedure  on  the  mass  and  momentum  equation  yields 


P2u2  sin  8  +  p2v2  cos  8  =  cf 


P2u2  sin  8  +  p2v2  cos  8  u2  +  P2  sin  8  =  c2 


p2u2  sin  6  v2  +  p2v2  cos  8  v2  +  p2  cos  6  =  c3 


(36) 

(37) 

(38) 


and  Equation  35  may  be  written  as 


P2u2  sin  8  h2  +  p2v2  cos  8  h2  +  p2u2  sin  8 


u2^ 


+  p2v2  cos  8  +  P2u2  sin  8  +  p2v2  003  & 


V2 


v22  = 


-  C/ 


(39) 


where  the  c  values  are  known  from  upstream  properties.  More  simply 
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p2u2  sin  B  +  P2v2  cos  B  =  ci 


(40) 


CjU2  +  P2  sin  3  =  C2 


(41) 


c|v2  +  ^2  cos  B  =  c3 


(42) 


and 


h2  +  -j-.  +-T-=  — 


u2 


c4 


(43) 


This  is  a  set  of  four  non- linear  algebraic  equations  which  were  solved 
for  the  downstream  properties  u2,  v2,  p2*  and  h2,  the  density  being 
determined  from  the  thermal  equation  of  state. 


IV.  DISSIPATION  MODEL 


Energy  loss  from  a  high  temperature  jet  is  chiefly  through  the 
mechanisms  of  conduction  and  radiation.  The  dissipation  term  required 
for  the  right  hand  side  of  Equation  20  Is  one  In  which  no  derivatives  of 
the  dependent  variables  appear.  The  usual  conduction  term,  therefore, 
cannot  be  utilized.  However,  there  are  several  radiation  models  that 
fulfill  this  requirement. 

in  general,  radiation  from  high  temperature  plasmas  is  con¬ 
sidered  to  be  the  continuum  due  to  free-free  and  free-bound  transitions 
of  the  electrons.  However,  free-bound  radiation  occurs  only  where  there 
is  recombination,  which  occurs  when  the  flow  is  near  equilibrium. 

At  pressures  In  the  range  of  0,01  and  0.5  atmospheres,  it  is 
not  uncommon  to  consider  the  flow  as  frozen  to  recombination.  Even  if 
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the  flow  Is  nearer  to  equilibrium  at  the  higher  pressures,  It  reduces 
computation  considerably  to  assume  frozen  flow.  This  allows  the 
investigator  to  study  trends  and  the  effects  of  various  parameters  on 
the  flow  with  but  a  fraction  of  the  computer  time. 

Assuming  the  flow  frozen  to  recombination,  the  radiation  would 
be  due  to  free-free  interactions.  Reference  (3)  gives  a  relation  due 
to  Spitzer  for  a  fuiiy  ionized  argon  gas  with  free-free  radiation  only  as 

Q  =  -  1.569  X  l(T34  n_  £  n.Cz.)2!1*  w-cuT3  (44) 

e  I  i  1 

However,  this  relation  gives  smaller  values  than  those 
measured  by  References  (4)  and  (5)  .  Multiplying  this  relation  by  15 
gives  a  better  approximation  to  these  results. 

Regardless  of  whether  this  approximation  Is  good  or  bad,  and 
there  is  enough  confusion  in  the  literature  to  suggest  that  It  is  as 
good  as  any.  It  serves  to  Illustrate  the  use  of  the  solution. 

V.  GAS  MODEL 

Following  Reference  (6),  for  an  ionized  gas  consider  the 

reaction 

Ar  ^  Ap+j  t  e  (45) 

where  Ar  is  an  ionized  particle  and  Ap+j  Is  the  next  higher  Ionized  par¬ 
ticle.  The  equilibrium  constant  for  the  reaction  can  be  expressed  as 
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nr+l  "e 


1 2ir  mekT  \ 

I  h2  I 


2  zrH 


(46) 


Here  Zr  is  the  electronic  partition  function  corresponding  to  particle 
Ar,  and  Zr+|  Is  the  partition  function  for  the  next  higher  ionization. 
The  parameter  ©r  is  the  characteristic  temperature  for  Ionization.  As 
the  gas  Is  assumed  to  be  a  mixture  of  perfect  gases,  the  equation  of 
state  may  be  written  for  three  levels  of  ionization  as 


=  no  +  n|  +  n2  +  n3  +  n0 
and  the  conservation  of  charge  is  stated  as 


(47) 


ne  =  n,  +  2n2  +  3n3  (48) 

Equations  46,  47,  and  48  are  a  system  of  five  non-linear 
equations  with  five  unknowns  which  were  solved  for  the  particle  densities. 
Thus  given  P  and  T,  density  is  returned  as 


p  =  mA  <no  +  +  n2  +  n^)  (49) 

The  energy  of  a  particle  is  divided  among  its  translational, 
rotational,  vibrational,  and  electronic  modes  (see  Reference  (7)).  This 
energy  depends  on  the  partition  function  for  each  mode.  In  general. 


3UnQr) 


Wf ~  =  kj2  '  "IT  t  Ii 


(50) 


where  is  the  total  energy  for  a  particular  species,  Ir  is  the 
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Ionization  energy  per  particle,  and  Np  Is  the  total  number  of  particles 
present,  and 


«  *  ^trans.  *  «rot  "  ®vtb  *  z  (5,) 

For  a  monatomic  gas,  the  rotational  and  vibrational  modes  are  not  present. 
Therefore,  Qro+  =  Qvjb  =  I,  and 


3 (An  Qfrans. ^  3 (in 

- §Y - +  -yr 


2r> 


t  I. 


From  statistical  mechanics. 


and 


^trans. 


2TrmkT  3/2 


-9|/T  -82/t 

zr  ■  CgD  +  g|  e  +  92  e  +  .  ..3r 


where  0j,  02  ...  ,  are  the  characteristic  temperatures  for  electronic 
excitation.  The  symbols  gQ,  gj,  ...  etc.,  represent  the  degeneracy 
factors  of  the  energy  levels  of  the  atomic  structure.  Specifically  for  a 
monatomic  gas 


gj  -  2J,  t  i 


where  Jj  are  the  inner  quantum  numbers  as  given  in  Reference  (8). 
Equation  50  then  becomes 


Er  3 


kT  + 


r  -0j/T  -02/T 

k[g|8|e  +  g282  e  + 

[g0  +  9|e'el/T  +  g2e_02/T  + 


t  I 


r 


(52) 
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Since  for  a  particular  species 


prV  =  NrkT 


and  the  enthalpy  Is  defined  as 


Hr  =  Er  +  prV 


The  ’’enthalpy  per  particle”  becomes 


Hr  5 
TC  =  7  kT  + 


-e./T 


-02/T 


k  [  9| 9|  e  +  92e2  6  + 

- 71 - 


■_v3.r  +  i. 


The  contribution  of  this  specie  to  the  specific  enthalpy  of  the  mixture  is 

Nrf5  k£V.I  3r.le-9r,«/T  ■ 

hr  =  BjT  kT  + - Tr - +  Ir_ 

where  Mj  is  the  total  mass  of  the  mixture.  By  summing  over  the  several 
species,  the  enthalpy  of  the  mixture  becomes 


h  = 


I 


i  t[i 


kT  +  kT 


ne  ^  kT 


(53) 


where 


Z  1  =1=0  for  electrons. 
©  © 


VI.  MACH  DISC 


The  location  of  the  so  called  Mach  disc  is  Important  to  many 
observers  as  its  attendant  subsonic  region  cannot  be  analyzed  by  the 
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characteristics  method.  Also,  the  method  of  characteristics  cannot  pre¬ 
dict  the  location  of  the  disc,  I.e.  the  location  of  the  embedded  oblique 
shock  Is  determined  until  It  Intersects  the  axis.  The  Mach  disc  occurs 
at  some  point  before  the  Intersection. 

Several  methods  have  been  offered  to  determine  the  Mach  disc 
location.  A  review  of  these  methods  Is  presented  In  Reference  (9).  A 
method  suggested  by  Eastman  and  Radtke(IO),  Is  simple  and  according  to 
Reference  (9)  is  reasonably  accurate  for  plume  exhausting  Into  still  air. 
They  hypothesize  that  the  Mach  disc  occurs  at  the  point  where  the  static 
pressure  behind  the  oblique  shock  reaches  a  minimum.  They  give  no  reason 
other  than  observation  to  support  their  contention.  Nevertheless,  their 
predictions  fall  relatively  close  to  other  more  sophisticated  predictions. 
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CHAPTER  111 

NUMERICAL  PROCEDURES 
I.  INTRODUCTION 


in  this  chapter  the  method  of  characteristics  will  be  formu¬ 
lated  for  machine  computation.  The  choice  of  procedure  is  largely 
subjective,  and  the  efficiency  of  one  procedure  over  another  is  not 
readily  apparent. 


II.  GENERAL  CONSIDERATIONS 

As  the  solution  of  the  characteristic  equations  is  an  initial 
value  problem,  the  flow  properties  along  some  starting  line  must  be 
known.  That  Is,  at  each  point  u,  v,  p,  and  h  must  be  specified. 
Arbitrarily,  the  starting  line  Is  chosen  to  be  the  exit  plane  of  the 
nozzle. 

The  computation  begins  at  the  starting  line  point  on  the  nozzle 
lip  as  In  Figure  3.  The  nozzle  exit  pressure  and  the  ambient  pressure 
exist  at  this  point.  Using  pressure  as  an  independent  variable,  the 
Prandt i -Meyer  relations  are  used  to  compute  a  set  of  properties  for  the 
point  for  each  arbitrary  increment  in  pressure  from  the  exit  pressure  to 
the  ambient  pressure.  The  next  point  down  on  the  starting  line  is  used 
in  conjunction  with  the  nozzle  lip  point  to  compute  a  third  point  down¬ 
stream.  The  nozzle  lip  point  is  treated  as  though  its  properties  were 
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Boundary—^ 


Fig.  3  Typical  Expansion  Network 
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those  existing  at  the  nozzle  exit. 

Using  this  new  point  and  the  nozzle  point  with  properties 
resulting  from  its  first  pressure  increment,  a  fourth  point  and  its 
attendant  properties  are  computed.  As  each  point  is  obtained,  its  pro¬ 
perties  are  stored  in  the  computer  for  further  use.  When  the  boundary 
Is  encountered,  i.e.  the  pressure  matches  the  boundary  pressure,  its 
properties  are  computed  and  stored.  This  storage  array  then  fully 
defines  a  complete  left-running  characteristic  curve.  Using  the  next 
point  on  the  starting  line  and  the  points  just  computed,  a  new  left¬ 
running  characteristic  is  computed,  and  so  on  until  the  starting  line 
point  on  the  axis  has  been  used.  The  next  characteristic  begins  from  an 
axis  point  downstream  from  the  starting  line.  It  is  computed  using  the 
method  of  reflection,  whereby  a  point  in  the  upper  half  plane  is  "re¬ 
flected"  In  the  lower  half  plane.  The  original  point  and  Its  "image"  are 
used  to  compute  the  point  on  the  axis. 

The  computation  continues  until  a  point  is  computed  that  is 
closer  to  the  beginning  of  a  left-running  characteristic  than  the  pre¬ 
vious  point.  This  indicates  that  a  "fold  back"  is  occurring.  At  this 
point  the  oblique  shock  wave  Is  considered  to  begin. 

As  in  Chapter  I,  the  field  above  the  point  of  foldback  is 
considered  the  upper  field  and  the  fiefd  below,  the  lower  field.  The 
field  properties  within  the  foldback  are  doubly  defined.  Guessing  a 
possible  shock  angle,  using  the  properties  of  the  point  as  defined  by 
the  lower  field,  and  substituting  Into  the  oblique  shock  relations,  a  set 
of  flow  properties  result  that  should  match  the  properties  of  the  point 
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as  defined  by  the  upper  field.  If  there  is  no  match,  another  shock  angle 
must  be  assumed,  and  so  on  until  a  matching  has  occurred. 

The  computation  then  proceeds  as  before  except  that  two  sets  of 
flow  properties  are  stored  for  each  shock  point.  The  computation  Is  seif 
terminating  when  the  shock  Intercepts  the  axis  or  when  the  shock  becomes 
so  strong  that  the  downstream  flow  becomes  subsonic. 

III.  TYPICAL  FIELD  POINT  CALCULATION 


In  Figure  4  the  position  and  flow  properties  of  points  I  and  2 
are  known  either  from  previous  computation  or  from  Initial  conditions. 

Putting  Equations  15,  18,  19,  20,  and  22  In  finite  difference 

form. 


^3 

“3 


- -  I—  2  -  2 

U | v |  +  a(yU|  -  8|  +  Vj 


u j2  $  a | ^ 


(54) 


(55) 


(56) 


4  -  X|  (57) 


P4  =X2  (58) 
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Point  to 


Fig.  4  Locating  of  Third  Point  from  Two  Known  Points 


AEDC-TR-69-184 


and 


P3V3V4  +  p3u3u4  +  p4  =  X3  (59) 

Vjv4  +  +  h4  =  X4  (60) 


X3  =  p3  +  P3V3V3  +  P3U3U3 

and 

X4  =  (I/^jUj)^  (Z4-Z3)  ±h,+  v3*3 
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Equations  55] and  56  are  used  to  find  the  tentative  position  of 
point  4.  Except  for  F,  the  barred  vaiues  are  initially  taken  to  be  the 
values  at  points  I  and  2.  The  F  values  are  averages  of  the  radial 
position  of  points  I  or  2  with  the  tentative  radial  position  of  point  4, 
i.e.  F|  =  h  <r|tr4).  The  streamline  through  points  3  and  4  Is  initially 
assumed  to  have  the  average  slope  of  the  streamlines  through  points  I 
and  2.  This  slope  and  the  position  of  point  4  initially  determines  the 
position  of  point  3.  The  properties  of  point  3  are  determined  by  linear 
interpolation  between  points  I  and  2.  Equations  57  to  60  along  with  the 
equation  of  state  constitute  a  set  of  linear  algebraic  equations  that 
are  solved  for  the  properties  of  point  4. 

The  barred  symbols  are  now  assigned  the  values  of  the  average 
of  the  properties  of  the  points  I  and  4,  2  and  4,  and  3  and  4,  i.e., 

IT |  =  %(U|  i-  u4) 

“3  ■  ^!u3  +  V 

“2  =  HlU2  +  V 

and  so  on. 

The  position  of  point  4  is  computed  again.  The  slope  of  the 
streamline  is  now  computed  from  Equation  54. 

The  position  of  point  3  Is  again  determined  In  the  same  manner, 
and  the  set  of  equations  is  solved  once  more.  This  iteration  Is  continued 
until  consecutive  values  of  the  properties  of  point  4  agree  within  some 
a  I lowable  error. 
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IV.  TYPICAL  PRANDTL-MEYER  EXPANSION 

Consider  Figure  5.  Multiple  sets  of  properties  are  assumed  to 
exist  at  point  2,  the  point  on  the  nozzle  lip.  The  pressure  is  allowed 
to  vary  continuously  from  its  value  at  the  nozzle  exit  to  the  ambient 
pressure.  The  expansion  fan  is  constructed  in  the  following  manner. 

Point  4  is  computed  as  a  typical  field  point  using  the  pro¬ 
perties  as  given  on  the  starting  line. 

Before  computing  point  4',  however,  the  Prandt l-Meyer  relations. 
Equations  27  and  28  along  with  Equation  26  must  be  used. 

The  pressure  differential  between  the  ambient  pressure  and  the 
nozzle  exit  pressure  Is  arbitrarily  broken  into  several  increments 
depending  on  network  size  desired.  Using  some  integrating  technique  such 
as  Runge-Kutta,  Equations  27  and  28  are  integrated  over  a  pressure 
increment,  determining  the  new  velocity  at  point  2.  Equation  26  yields 
the  enthalpy  directly.  Using  these  new  properties  for  point  2  and  the 
properties  of  point  4,  point  4’  is  readily  determined  as  a  typical  field 
point  calculation. 

This  procedure  is  continued  until  the  ambient  pressure  is 
reached  at  point  2.  The  boundary  point  is  then  computed  as  a  typical 
boundary  point. 


V.  TYPICAL  BOUNDARY  POINT  CALCULATION 

Figure  6  could  be  considered  as  "half"  of  Figure  4,  page  25. 
That  is  to  say,  since  the  boundary  is  formed  by  the  streamline,  the 
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Fig.  5  Typical  Expansion  Fan  Construction 
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Fig.  6  Locating  of  Second  Boundary  Point  from  a  Previously 
Known  Boundary  Point  and  Field  Point 
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characteristic  line  between  points  2  and  4  in  Figure  4,  page  25,  lies 
outside  the  region  of  interest.  Two  boundary  conditions  replace  the  two 
relations  defining  this  characteristic  line  and  its  properties. 

One  condition  is  simply  that  the  boundary  is  a  streamline. 

The  second  condition  Is  that  for  a  pressure  boundary,  the  pressure  is 
specified  at  all  points  along  the  boundary  streamline,  i.e. 

p  =  f(z,r,u,v,h) 

For  the  case  considered  here,  the  boundary  is  a  free  pressure  boundary, 

I  .e. 


p  =  constant 


The  required  difference  relations  are 


(61) 


(62) 


/ai-1 

/  Ar  \ 

LU| lAz#LI 

-  vl. 

UiV4  - 

.u'  Nli  '  V|. 

V|U4  =  Xg| 


(63) 


and 


p3v3v4  +  p3u3u4  =  f>,  v,  ^  -t- u,  Uj 

V3V4  +  U3U4  +  h4  =  XB3 


(64) 

(65) 
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where 


and 

XB3  =  X4 

The  expressions  for  X|  and  X?  are  unchanged  from  those  in  the  typical 
field  point  calculation. 

The  solution  procedure  is  similar  to  that  of  the  typical  field 
point,  except  that  point  3  Is  fixed.  Point  4  Is  tentatively  found  from 
Equations  61  and  62,  then  the  linear  set,  Equations  63  to  65,  are  solved 
for  the  properties  at  point  4.  New  barred  values  are  assigned,  and  so  on 
following  an  iterative  procedure. 

VI.  TYPICAL  POINT  ON  THE  AXIS 


This  being  an  axially  symmetric  problem,  the  calculation  of 
axial  points  Is  quite  simple  using  the  method  of  reflection.  From 
Figure  7,  point  I  is  reflected  In  the  lower  half  plane  as  point  2,  I.e. 
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Fig.  7  Locating  of  Paint  on  Axis  Using  Reflected  Point 
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P2  =  P| 
and 


The  calculation  procedure  is  then  that  of  a  typical  field  point. 

VII.  TYPICAL  SHOCK  POINT 

The  starting  of  the  Interior  shock  is  signaled  when  the  con¬ 
dition  depicted  in  Figure  8  occurs,  i.e.  when  point  4  falls  behind  point 
I,  meaning  that  two  right-running  characteristics  have  crossed.  Follow¬ 
ing  Reference  ( 1 1 ),  the  left-running  characteristic  curve  B  is  completed 
to  the  boundary  and  point  4  Is  deleted  from  the  field.  Characteristic 
curve  C  Is  then  completed  up  to  the  virtual  point  (point  VP).  Point  I 
Is  considered  as  the  first  shock  point. 

A  guess  of  the  shock  angle  8  is  made  such  that  the  tentative 
shock  line  fails  to  the  right  of  the  characteristic  line  between  points 
I  and  VP.  Guessing  8  such  that  the  shock  line  fell  to  the  left  results 
In  an  Incorrect  solution,  i.e.  no  solut-ion  of  the  oblique  shock  rela¬ 
tions  can  be  found  that  has  an  entropy  increase.  This  Is  because  the 
limits  of  a  shock  angle  are  between  that  of  a  Mach  line,  here  a  charac¬ 
teristic,  and  a  shock  normal  to  the  upstream  flow. 

The  intersection  of  the  shock  line  with  characteristic  C  is 
found  and  its  properties  obtained  by  interpolation.  These  properties 
are  the  upstream  properties  at  the  shock  point  (point  SPI). 
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Fig.  8  Matching  of  Upper  Field  to  Lower  Field  Through  A  Shock  Point 
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and 


The  shock  relations. 


P2u2  sin  6  +  p2v2  cos  9  =  P|U}  sin  9  +  p(V|  cos  6  =  Cj 


clu2  +  P2  sin  9  =  c | u |  +  pj  sin  9 
civ2  +  ^2  cos  ®  "  c|v|  +  Pj  cos  9 


l  u22  +  v22  i.  ,  ur  +  vr 
h2  +  2  »  h,  +  2 


ui2  +  Vj2 


are  then  solved  for  the  downstream  properties  at  point  SP2.  (Note  that 
points  SPI  and  SP2  have  the  same  position.)  The  constants  are,  of 
course,  determined  from  the  upstream  properties,  while  the  subscript 
(2)  here  refers  to  downstream  properties. 

For  this  to  be  a  solution,  these  properties  must  match  those 
predicted  by  the  upper  field  for  this  position. 

First,  construct  a  right-running  characteristic  through  point 
SP2  to  some  point  BP.  The  tentative  slope  is  given  by 


Ar  _ 
Az  = 


where  the  barred  quantities  are  Initially  assigned  the  values  of  the 
properties  at  point  SP2.  The  tentative  point  BP  Is  found  and  its  pro¬ 
perties  interpolated  for  In  characteristic  B.  The  barred  values  are 
then  assigned  the  value  of  the  average  properties  of  points  SP2  and 
BP.  A  new  slope  is  found  and  then  a  new  position,  and  so  on.  The 
Iteration  Is  stopped  when  succeeding  positions  of  point  BP  change  less 
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than  some  allowable  error.  The  barred  values  and  the  delta  values  are 
then  substituted  into  the  finite  difference  form  of  the  compatibility 
relation  for  a  right-running  characteristic,  i.e. 


If  the  substitutions  do  not  satisfy  this  equation,  another  guess  as  to 
the  shock  angle  must  be  made,  and  an  iteration  devised  such  that  this 
equation  is  finally  satisfied. 

After  convergence,  the  final  properties  of  point  SP2  are  used 
in  conjunction  with  the  properties  of  the  next  point  up  along  charac¬ 
teristic  B  to  compute  the  next  point  on  characteristic  C,  and  so  on  as 
before  until  the  boundary  is  reached. 
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CHAPTER  IV 
SAMPLE  PROBLEM 

I.  GENERAL  CONSIDERATIONS 


There  has  been  some  recent  experimental  work  done  in  the  area 
of  high  temperature  gasdynamlcs  using  an  arc  heater  as  an  energy  source 
and  argon  as  the  gas.  At  AEDC  work  is  progressing  on  a  study  to  determine 
the  feasibility  of  accelerating  a  high  temperature  jet  with  an  electric 
field  for  possible  wind  tunnel  application.  Though  the  primary  interest 
is  in  the  effects  of  the  magnetohydrodynamic  forces,  knowledge  of  the 
plume  structure  without  electrical  effects  is  of  considerable  importance 


II.  SOLUTION 

Workers  with  arc  heaters  have,  in  general,  been  interested  In 
jets  exhausting  from  a  sonic  nozzle  in  the  temperature  range  of  14,000 
to  20,000°K  at  a  pressure  of  the  order  of  0.01  to  I  atmospheres  into  an 
ambience  at  a  pressure  of  from  one-eighth  to  one-fiftieth  of  the  exit 
pressure. 

Figure  9  shows  the  plume  resulting  from  the  following  exit 
conditions  without  radiation  assuming  the  flow  is  frozen  at  the  exit 
equilibrium  conditions: 

i.  Mach  number  =  1.0 1 . 

.  2.  Exit  pressure  =  33.42  pounds  per  square  foot. 
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Starting  Conditions: 

Nozzle  Half-Angle  =  10  deg 
Exit  Pressure  =  33.42  lbf/ft2 
Ambient  Pressure  =  0.668  lb^/ft2 
Temperature  =  14,000°K 


4.75 


sn3.6' 

Shock 


Line  of  Reversed  Foldback 

—I _ I _ I 

0.6  0.7  0.8  • 


Fig.  9  Plume  with  Negligible  Radiation  Effect  Showing  Lines  of  Constant  Mach  number 
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3.  Ambient  pressure  =  0.668  pounds  per  square  foot, 

4.  Temperature  =  I4,000°K. 

5.  Nozzle  radius  =  0.02625  feet. 

Including  radiation  with  the  above  starting  conditions  had 
negligible  effect.  The  difference  was  so  small  that  It  was  Impossible 
to  plot. 

Consider  Figure  10,  the  Eastman  and  Radtke  method  for  locating 
the  Mach  disc,  as  discussed  In  Chapterll,  failed  in  this  case  since  the 
computer  program  terminated  before  the  minimum  point  on  the  pressure  plot 
occurred.  The  program  terminated  because  a  crossing  of  left  running 
characteristics  was  occurring  Indicating  that  another  shock  was  trying  to 
form  along  the  dashed  line  in  Figure  9.  The  program  cannot  handle  this 
reversed  foldback.  It  would  be  of  no  value  to  have  the  program  construct 
this  shock  as  in  the  actual  plume  the  Mach  disc  and  Its  attendant  sub¬ 
sonic  region  would  be  found  somewhere  In  this  region. 

In  order  to  find  a  significant  radiation  effect,  the  pres¬ 
sures  and  temperature  were  raised.  The  following  case  then  was 
run  both  with  and  without  radiation  to  demonstrate  the  non-lsoenergetlc 
solution: 

1 .  Mach  number  =  1.01. 

2.  Exit  pressure  =  2116.22  pounds  per  square  foot. 

3.  Ambient  pressure  =  264.53  pounds  per  square  foot. 

4.  Temperature  =  20,000°K. 

5.  Nozzle  radius  =  0.02625  feet. 

Here  again  the  flow  was  assumed  frozen  at  equilibrium  exit 
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conditions.  Though  here  it  is  not  possible  to  justify  frozen  flow  on 
physical  grounds,  it  becomes  Justifiable  when  the  primary  interest  is  in 
changes  in  the  major  features  of  the  flow  field  and  computer  time  is  at  a 
premium.  Thus,  though  the  program  can  handle  equilibrium  flow  with 
facility,  the  required  computer  time  differs  from  that  for  frozen  flow 
by  an  order  of  magnitude. 

Comparing  Figure  II  with  Figure  12,  the  radiation  effects  are 
apparent.  Though  the  plume  shape  and  the  location  of  the  embedded  shock 
change  little,  there  is  a  definite  shift  in  the  constant  Mach  number 
lines.  The  shift  Is  most  apparent  In  the  higher  density  region  of  the 
upper  field.  The  flow  near  the  axis,  however,  is  more  easily  analyzed. 
There,  as  would  be  expected,  the  energy  loss  due  to  the  radiation  raises 
the  Mach  number  at  a  particular  location.  This  is  so  since  the  radiation 
tends  to  lower  the  temperature  thus  lowering  the  acoustic  speed  while 
affecting  the  velocity  very  little.  The  upper  field  Is  more  complex  in 
that  a  streamline  will  pass  from  a  lower  number  to  a  higher  one  then  back 
to  a  lower  Mach  number.  In  general  it  can  be  said  that  when  passing  to 
the  higher  Mach  number,  the  flow  reaches  the  higher  Mach  number  sooner 
with  radiation.  Similarly  when  passing  to  a  lower  Mach  number,  the 
reaching  of  the  lower  Mach  number  is  delayed  somewhat  with  radiation. 

In  Figures  13  and  14,  the  Eastman  and  Radtke  criterion  was 
used  to  locate  the  Mach  disc.  The  method  could  not  be  applied  with 
assurance  in  Figure  13  as  the  curve  did  not  have  a  unique  minimum  point. 
The  curve  with  radiation  did,  however,  have  a  unique  minimum.  The 
apparent  effect  of  the  radiation  was  the  shifting  of  the  Mach  disc 
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Starting  Conditions: 

Nozzle  Half-Angle  «  10  deg 
Mach  Number  ■  1.01 
Exit  Pressure  -2116. 22  Ibftft2 
Ambient  Pressure  -  264.5  lbf/ft2 
Temperature  -  20, 000°K 


2.35 


z,  ft 


Fig.  11  Plume  without  Radiation  Showing  Lines  of  Constant  Mach  Number 
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Location  of 
Mach  Disc 


I  I  1 _ I - 1 _ I _ I _ U _ I _ I 

0  0.02  0.04  0.06  0.08  0.10  0.12  0.14  0.16  0.18  0.20 

z,  ft 


Fig.  13  Pressure  Distribution  Behind  Embedded  Shock  for  Case  with  No  Radiation 
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.04  0.06  0.08  0.10  0.12  0.14  0.16  0.18  0.20 

z,  ft 


Fig.  14  Pressure  Distribution  Behind  Embedded  Shock  for  Case  with  Radiation 
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downstream.  The  reason  for  this  effect  is  not  obvious. 

As  to  the  accuracy  of  the  method,  it  can  only  be  said  that 
the  results  with  no  radiation  agree  very  well  with  the  results  in 
Reference  (II)  and  other  method  of  characteristics  solutions  which  are 
too  numerous  to  list. 

Due  to  the  confusion  in  the  literature  as  to  just  what 
constitutes  a  good  radiation  model,  it  can  only  be  concluded  that  the 
results  agree  qualitatively  with  what  one  expects  to  occur.  That  is, 
the  direction  of  shift  of  the  constant  Mach  number  lines  and  the  increase 
of  radiation  effect  with  temperature  and  pressure  are  as  what  would  be 
expected. 
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CHAPTER  V 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  following  conclusions  are  made  as  a  result  of  this  study: 

1.  Even  though  the  inclusion  of  the  non-isoenergetlc  term 
in  the  energy  equation  and  the  inclusion  of  the  embedded 
shock  Into  the  method  of  characteristics  solution 
greatly  complicates  the  numerical  procedures  Involved, 
the  resulting  computer  program  is  much  more  flexible  and 
provides  a  more  physically  correct  solution. 

2.  The  use  of  the  non-lsoenergetlc  solution  has  been 
demonstrated  by  comparison  of  flows  with  and  without 
radiation.  These  comparisons  show  that  radiation  has 
little  effect  on  the  flow  except  at  relatively  high 
density,  i.e.  in  the  temperature  range  15,000  to 
20,000°K,  the  radiation  is  not  important  except  at 
pressures  of  the  order  of  one-half  atmosphere  or  above. 

The  following  recommendations  are  made  for  future  studies: 

1.  A  study  should  be  made  to  correlate  the  non-isoenergetlc 
solution  with  experiments  leading,  perhaps,  to  a  better 
dissipation  model . 

2.  The  computer  program  should  be  further  developed  so  as  to  handle 
walls  and  center  bodies,  thus  making  the  program  useful  In 
supplementing  nozzle  studies  where  real  gas  effects  are  Important. 
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